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ABSTRACT 



We have performed MHD simulations of Alfven wave propagation along an 
open flux tube in the solar atmosphere. In our numerical model, Alfven waves 
are generated by the photospheric granular motion. As the wave generator, we 
used a derived temporal spectrum of the photospheric granular motion from G- 
band movies of Hinode/SOT. It is shown that the total energy flux at the corona 
becomes larger and the transition region height becomes higher in the case when 
we use the observed spectrum rather than white/pink noise spectrum as the 
wave generator. This difference can be explained by the Alfven wave resonance 
between the photosphere and the transition region. After performing Fourier 
analysis on our numerical results, we have found that the region between the 
photosphere and the transition region becomes an Alfven wave resonant cavity. 
We have confirmed that there are at least three resonant frequencies, 1, 3 and 5 
mHz, in our numerical model. Alfven wave resonance is one of the most effective 
mechanisms to explain the dynamics of the spicules and the sufficient energy flux 
to heat the corona. 

Subject headings: Sun: corona — Sun: granulation — Sun: oscillations — Sun: 
photosphere — Sun: solar wind — Sun: transition region 
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Introduction 



From the time the solar corona was found to be at a temperature of more than 1 



million K ( lEdlen 



19431 ). the search for the process that heats the corona, so called coronal 
heating problem, has become one of the most important problems in the solar physics. 
It is known that mechanical energy transport is required to heat the corona because the 
photospheric temperature is at most 6,000 K. Since large amount of mechanical energy flux 
can be generated by the photospheric convection motion, the essential problem is how to 
transport the mechanical energy to the corona and how to dissipate it in the corona. The 
mechanisms that heat the actual corona are yet to be elucidated, since it is difficult to 
self-consistently include all the energy transport process. 

It is required that the heating rate of any coronal heating mechanisms must balance 
the energy loss by radiative cooling and thermal conduction in order to maintain the 1 
milli on K corona. The energy flux necessary to heat the quiet corona is 3 x 10^ erg cm~^ 



s ^ ( IWithbroe fc Noyes 



19771 ). There are two mechanis ms that are widely believed to heat 



the c orona. There are the Alfven wa ye heating ra odel ( Uchida &: Kaburak 



19741 ). and the the nanofiare model ((Parker 



1974; 



Wentzel 



19881 ). Both models explain how to release 



Poynting flux generated by the photospheric convection motion. 

In the Alfven wave heating model, magnetoacoustic waves driven by the convection 
motion propagate upward to produce shocks. Damping length of compressional waves such 
as fast or slow mode waves is very short because of the highly stratified atmosphere in the 
Sun. Moreover the diffraction due to the gradient of Alfven speed prevents compressional 
waves from propagating to the corona. Therefo re it is consid ered that it is not possible to 



19811 ). On the other hand, Alfven 



heat the corona by these compressional waves (iHoUwegi 
waves are non- compressional waves, so they can reach the corona even if the atmosphere 
is highly stratified. Damping length of Alfven modes is defined by various dissipation 
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Suzuki fc Inutsuka 



2005 



( iHeyvaerts fc Priest 



(HoUwee et al. 


198^ 


)■ 


Ku 


doh & Shibata 


1999 


resonance ( 


lonson 


1978 


), and phase mixing 



19831 ). Alfven wave resonance also occurs when a flux tube between 



the photosphere and the chro mosphere acts as an Alfve n wave resonant cavity, this can 



happen even if it is not closed ( jHollweg 



1981 



Leroy 



19811 ) ■ Since a large amount of the wave 



flux is reflected by the transition region, standing waves appear in the region between the 
photosphere and the transition region. These standing waves can be amplified depending 
on the wave frequency excited at the wave source and they can eventually be dissipated 
nonlinearly. It is not understood yet which processes can act efficiently in the real solar 
corona. 

Various studies have been ma de to investigate the n onlin ear propagation of A l fven 



waves using numerical simulations. 



Moriyasu et al. 



fl2004h and Morivasu fc Shibatal fl2004h 



studied the dissipation process of Alfven waves in coronal loops. In their models, shear 
Alfven waves are driven by footpoint motion, and are coupled nonlinealy with fast and slow 
mode waves during propagation in the chromosphere and the corona. They suggested that 
Alfven waves can create 1 million K corona w hen the r.m.s veloc i ty of the footpoint motion 



is greater than 1 km s ^. For open flux tubes. 



Kudoh fc Shibatal (119991 ) showed that Alfven 



wave heating model can also explain the sp icule motion as wel 



energy flux to heat the corona. Moreover, 



Suzuki fc Inutsuka 



as it can tra nsport sufficient 



(12005 . 



20061 ) showed that 



nonlinear low-frequency Alfven waves can simultaneously heat the corona and drive the fast 
solar wind. 

One issue with the above simulations is the wave generation in the photosphere. 
Because ther e were only few obs ervations on the time series of transverse motion at the 



photosphere (ITarbell et al 



19901), the previous simulations assumed an artificial velocity 



spectrum for wave generation. As both the total power and the shape of the velocity 
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spectrum can affect the wave generation and the heating process, it is necessary to use 
observational data to justify whether Alfven waves are able to heat the corona. 

It is possible to observe turbulent convection motion in the photosphere by Doppler 
analysis. Temporal velo city profiles for the vertical motion have bee n studied for the 



sake of helioseismology ( jHarvey et al. 



19931 ). ICarlsson fc SteinI ( 119971 ) used the observed 
Doppler velocity as a generation mechanism of acoustic waves and reproduced the Ca II 
spectrum to compare with observations. On the other hand, generation of Alfven waves 
requires transverse motion in the photosphe r e. Local correlatio n tracking (LCT) method 



( iNovember fc Simon 



1988 



Title et al 



1989 



Berger et al. 



19981 ) reveals the horizontal 



motion of the granules. T hough many studies have been done on the s patial distribution 



of the horizontal velocity ( iNovember fc Simon 



1988 



Berger et al 



19981 ). few studies have 



1990l ). This is because the 



been performed on the temporal velocity profiles ( iTarbell et al. 
LCT method is greatly affected by the atmospheric seeing. However, a seeing free data set 



can be obtained by Solar Optical Telescope, SOT ( iTsuneta et al 



( iKosugi et al. 



20081 ) on board Hinode 



20071), so we can clearly understand the properties of the temporal profiles of 



the horizontal motion through LCT techniques. 

An important physical process in the Alfven wave dynamics is the spicule dynam ics. 
One of the mechanism of spicule formation is s. 



Shibata fc Suematsu 



Suematsu 



1982 



Suematsu et al. 



1982 



ow mode shock moc 


lei ( 


HoUwee 


1982 




Sterline; & HoUwee 


1988; 


Sterling 


2000) 



fll985l ) applied the slow shock model to fibrils and threads in a loop g eometry. 



Leak a ge of p-mode from the photosphere even 



1990 



De Pontieu et al. 



2004 



tually produce the slow mode shock (ISuematsu 



Hansteen et al. 



20061 ) . These shocks lift up the transition 



region resulting in spicule motions. However, spi cule height cann ot reach the observed value 



(5,0 00-7,000 km ) 



and 



Saito et al. 



only by the leakage of p-mode ([Suematsu 



1990) 



Kudoh fc Shibatal mm 



(l200l[ ) considered the slow mode generation through nonlinear Alfven wave 
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mode conversion. Thanks to Ca II n iovies observed by Hinod e satellite, detailed dynamics o: 



the s pi cules have b een u nders tood ( 



2009h . 



KimetaL 



(120081 ) and 



De 



He et al 



'ontieu et al. 



2007a 



Rouppe van der Voort et al. 



(|2009[ ) estimated the minimum wavelength of kink 



Alfven waves from horizontal motions of spicule. These observation can impose many 
restrictions on coronal heating models. 

The purpose of the present study is to include a realistic horizontal velocity spectrum 
in the magnetohydrodynamic simulations and justify coronal heating by Alfven waves in 
the solar atmosphere. In section 2, we will describe observation method to measure the 
horizontal photospheric motion. The basic assumptions and our numerical model will be 
explained in section 3. In section 4 we describe our numerical results and discussion will be 
included in section 5. Finally, we will summarize our conclusion in section 6. 



2. Observations & Data Reductions 



One of the main purposes of the Hinode satellite is to solve the coronal heating 



problem. Among the three telescopes 



( iTsuneta et al. 



2008 



Suematsu et al. 



equipped on Hinode, Sola r 



2008 



Ichimoto et al 



2008 



Optical Telescope 



Shimizu et al. 



(SOT) 



20081) is an 



optical telescope which has high spatial resolution (0."2-0."3) depending on wave length. 
Since observations from space are not affected by atmospheric seeing, Hinode makes it easy 
to obtain a seeing free data set over a long time span. 

In this study, we used fourteen data sets from quiet regions observed with G-band filter 
between 2006 October 31 and 2007 December 29. We described the properties of each data 
set in Table [H The mean time cadence of each data set is less than 35 second. Basically, we 
selected the data sets whose duration is longer than 70 min. Since every data set is selected 
to locate within ±100 arcsec from the disc center, where the inclination angle is less than 
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6 degrees. Pixel resolution of the CCD is 0/'054 and FOV is larger than 512 pixel x 512 
pixel. Figure [T] shows the typical G-band image used in our analysis. 

Dark current subtraction and flat-fielding were applied in the standard manner for 
all of the G-band images. There are overall displacements between two images due 
to the solar rotation, tracking error, and satellite jitter. We have calculated the cross 
correlation between two consecutive images and derived the displacement that had the 
highest correlation coefficient. Using the displacement, we have adjusted all the data to an 
accuracy of sub pixel scale. 

Generally, C-band images also have intensity fluctuations from p mode as well as 
granular motion. In order to reduc e the contribution from p or f modes, we have applied so 



called sub-sonic flltering technique (ITitle et al. 



19891 ) to our data set. The sub-sonic flltering 



technique removes the power where u/k > Vph in Fourier space so that the information 
moving faster than Vph does not appear in the movie. In the present study, we have selected 

ph 



Vnh = 5 km s ^ 



In order to obtain the photospheric horizontal velocity, we have 



applied the LCT 



Berger et al. 



(11998h 



technique to the C-band movies. The LCT techniques developed by 
uses rectangular subflelds or "tiles" as a tracer and determines the displacement of tiles 
between two images. Since the horizontal velocity derived by the LCT technique usually 
depends on the size of the tracer, we have flxed the tile size to be 0."4. Using the Fast 
Fourier Transform technique, we have estimated the one-sided power spectrum density, P^, , 
from the time series of the photospheric horizontal velocity. The one-sided power spectrum 
density is deflned as, 

<V?>= I P„dv. il] 



where V±_ is the horizontal velocity, and z/^m, i^max are the minimum and maximum 
frequencies constrained by the observational sampling time and duration. <> means the 
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temporal average. Power spectrum densities were first calculated at each fixed point, then 
we averaged over the FOV. Figure [2] shows an example of the power spectrum density in 
the region # 1 listed in Table [TJ The solid line represents the mean power and the dotted 
lines above and below the solid line represent the range of 1 sigma. 

We have also used Ca II movie observed by Hinode/SOT on 22-Nov-2006 in order to 
compare apparent motion of spicules with our numerical simulation. The target region is 
the west limb over the quiet region. Pixel resolution of each image is 0."054 and cadence is 
8 sec. 



3. Numerical settings 



Kudoh fc Shibatal f ll999h . 



The basic concepts of our simulations are the same as those of 
We consider the nonlinear propagation of Alfven waves driven by photospheric convection. 
The whole system is described by MHD equations with polytropic index 7 = 5/3, including 
gravity and empirical chromospheric cooling. All physical variables are supposed to be the 
function of time t and the distance s along the flux tube measured from the photosphere. 
Vector variables are composed of the component along the flux tube (poloidal component) 
and the toroidal direction. We suppose the mean molecular weight to be 1, independent of 
the plasma status. Moreover the poloidal configuration of the flux tube is assumed to be 
rigid, though the toroidal components of velocity (V^) and magnetic field (B^) are allowed 
to change with time. Above assumptions lead the following basic equations. 



(2) 



PVS^ + Pg 



Bl 



in 



A 



pV^Vs 



B^Bs 
An 



^3/2) =ApL 



trqi 



(4) 



di 



(yAB^^ + — ([B^Vs - BsV^] = 0, 



(5) 



d_ 

di 
d_ 

ds 



P|V| 



+ 



7-1 



+ 



|B| 



57r 



7 



47r 



B^BsVff, 



An 



A 



dz 



-LradA - pVsQ—A + pV^^^ALtrq, 



(6) 



p 



pksT 



rrir. 



(7) 



where r is the radius of the flux tube, A is the cross section of the flux tube, p is the 
mass density, Pg is the gas pressure, T is the temperature, Vg is the poloidal component of 
the velocity, is the toroidal component of the velocity, L^rq is the photospheric torque 
used as a wave driver, Bg is the poloidal component of the magnetic field, and B^ is the 
toroidal com ponent of the ma gnetic field. We include the chromospheric cooling Lrad 

(119971 1 only when the plasma temperature is less than 8 x 10^ K. 



according to 



Hori et al. 



This method involves two processes. For plasmas whose tempera ture is less than 4 x 10 



we implement 



Sterling et al. 



the e mpirical cooling for optically thick plasmas ([Anderson Sz Athay 



K, 



1989 



19931 ) 
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L,„rf = 4.9 X 10V(s,t), ergcm-^K-i. (8) 
For plasmas whose temperature is more than 4 X 10"^ K and less than 8 x 10^ K, we use the 



cooling for optically thin plasmas ( iHildner 



1974 



Sterling et al. 



19931 ) 



Although flux concentrations at the photosphere are on the order of the pressure scale 
height, they rapidly expand in the upper atmosphere due to exponential decrease of gas 
pressure. If we assume mechanical equilibrium between the magnetic pressure of the flux 
tube and the ambient gas pressure, plasma beta is nearly 1 in the photosphere and the 
chromosphere. In o rder to reproduce such a situation, we assume the same structure of the 



flux tube as that of 



Kudoh k Shibatal fll999h . 



where 



cos ads, (9) 
sinarfs, (10) 

a = at + {ar- at)fn, (11) 

at = — arctan [{z/zd) ] , (13) 

/„ = -^{tanh[(^ - 0.2^rf)/(0.1^,)] - 1}. (14) 

Hq is pressure scale height of the photosphere (150 km), Zd is equal to 12 Hq (1800 km) 
, and r is the radius of the flux tube {A oc r^). The poloidal fleld Bg is determined by 
Bs = BQlro/r)"^, where is the radius of the flux tube at the photosphere {tq = Hq = 150 
km), and Bq is deflned so that plasma beta equals to 1 in the photosphere (Fig. E]). 

The initial atmosphere is supposed to be magnetohydrostatic equilibrium with constant 
gravity {g = 2.7 x 10^ cm s~^). The temperature is taken to be as follows. 
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- = l + -(a.-l)(l+tanh^^j, (15) 

where Oc = 300, and Tq is the photospheric temperature (5000 K). The density at the 
photosphere is supposed to be 2.53 x 10^'' g cm~^. Figure S] shows the profile of the plasma 
beta with respect to height in our model atmosphere. 

The foot point of the fiux tube is subjected to an external torque [Ltrq) so as to drive 
the Alfven waves. The formulation of the external torque is 



L., = -xFoX^l-tanh^^^^J. (16) 

The amplitude of the torque, Fq, is determined so that the resulting toroidal velocity of 
the foot point becomes the given velocity time series that reproduces the given velocity 
spectrum. 

The numerical domain ranges from the photosphere to more than 900 Hq (= 1.35 xlO^ 
km). The grid size below 110 Hq from the photosphere is 15 km while the grid size over 
110 Hq becomes 1.01 times larger per grid. Since wave amplitude increases exponentially 
with height in the stratified atmosphere, we set the gravity to be above 900 Hq in order 
to reduce the effect of reflection at the upper boundary. 

Our numerical schem e is an HLLD approximate Riemann solver developed by 



Miyoshi &: Kusand ( 120051 ). Since the HLLD scheme has positivity in some situations, this 
scheme is more robust than a Roe-type linear Riemann solver. Moreover this HLLD scheme 
has good resolution and numerical efficiency. The spatial accuracy obtains 2nd order by 
using the TVD-MUSCL method while the temporal accuracy obtains 2nd order by using 
the Runge-Kutta method. We expect that a conservative scheme like t he HLLD scheme 



can ca pture the shocks better than the non-conservative scheme used in 



Kudoh fc Shibata 



( 119991 ). Conservative schemes are suitable for the situation when Alfven waves are converted 



into compressional waves resulting in a lot of shock formation. 



4. Results 

We have investigated the relation between the initial Alfven wave spectrum and the 
nonlinear behavior of the upper atmosphere. After assuming the model atmosphere, Alfven 
waves are generated by torque at the foot point of a magnetic flux tube. First of all, 
monochromatic waves that have 1 km s~^ amplitude at the photosphere are generated. 
Then, we have driven Alfven waves that have continuum spectrum with different wave 
amplitudes. Finally, the observed velocity spectrum is included in our numerical simulation 
as a wave driver. 

Figure [5] shows the dynamics of Alfven waves in typical case of our numerical simulation. 
Plots in each panel are stacked with time increasing upward in uniform increments of 8.0 s. 
Although we only generated toroidal motion, poloidal motion can be seen in the figure. This 
is the result of mode conversion from Alfven waves to slow mode waves. The slow mode 



wave s push up the transition region, which can be observed as spicules (IKudoh fc Shibata 
1999k 



4.1. Monochromatic wave 

We generated monochromatic Alfven waves, changing the wave period from 10 sec 
to 750 sec. For each period, the amplitude is fixed to be 1 km s~^ and we measured the 
resulting energy flux in the corona (15,000 km) and the height of the transition region. 

Figure [6] shows the relationship between the mean total energy flux measured at the 
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corona and the wave period. Total energy flux is described as follows; 



Ip{V^ + V,^) + -^^P, + P9z+^ 



B^Bs 
An 



V^. (17) 



The dotte d line represents the suffic ient energy flux to maintain the hot corona , 3 x 10^ erg 



cm ^ s ^ (IWithbroe fc Noyes 



19771 ). The vertical solid lines at each data point represent 
the error bar for 1 sigma. The energy flux starts to decrease rapidly when the wave 
period becomes larger than 500 sec. The time span, 500 sec, roughly corresponds to the 
Alfven wave travel time back and forth between the photosphere and the transition region. 
Therefore the waves whose period are longer than 500 sec are canceled out when they come 
back to the photosphere. The waves whose period are shorter than 60 sec also start to 
decrease their energy flux at the coronal height. In the case of the shorter period Alfven 
waves, the shorter period longitudinal waves are frequently produced through nonlinear 
mode conversion. Because such waves are large in number and have shorter damping 
length, a large amount of wave power is dissipated in the chromosphere. Wave periods that 
efficiently heat the corona range from 60 to 500 sec. One of the reason for that is the Alfven 
wave resonance between the photosphere and the chromosphere, which we discuss in detail 
in the following section. The energy flux profile with respect to wave period (Fig [6]) clearly 
shows resonant peaks at 300 and 450 sec. 

Figure [7] shows the relationship between the mean transition region height and the 
wave period. In this study the transition region is defined as the region of the atmosphere 
at a temperature of 8 x 10^ K. Similar to the mean energy flux, the waves whose period is 
shorter than 60 sec or longer than 500 sec cannot create the tall spicules that are frequently 
observed in the real Sun. The transition region height profiles with respect to wave period 
(Fig. [7]) also represent the resonant peaks at the 130, 300, 450 sec. 
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4.2. The waves that have continuum spectrum 

We have investigated the dependence of nonUnear Alfven wave propagation on the 
spectral shape of the initial waves. Two different wave spectrums, white noise and pink 
noise, are considered as wave generators. The power spectrum densities P^, defined by 
equation ([1]) of the white noise and pink noise are proportional to z/° and z/~^ respectively. 
The frequency range is li mited from 6 x 10"'^ to 5 x 10^^ Hz in the same way as that of 



Suzuki fc Inutsukal (120061 ). For each spectrum, the amplitude of the waves changes between 



0.1 and 5 km s ^. 

Figure [8] shows the dependence of the total energy fiux in the corona on the amplitude 
of the wave. The horizontal dashed line represents the sufficient energy fiux to heat the 
corona. The diamond and asterisk symbols show the case for white noise and pink noise 
respectively. Generally energy fiux at the corona increases with the wave amplitude. In 
white noise case, we need a wave amplitude more than 2 km s~^ to maintain the hot 
corona while in the pink noise case, the necessary wave amplitude is more than 1.3 km 
s~^. Pink noise waves heat the corona more efficiently than white noise waves because the 
pink noise waves include more power in the lower periods around 300 sec. Around these 
lower period ranges, Alfven wave resonant periods exist. Therefore large amount of energy 
fiux can propagate to the corona without refiection. In the case of pink noise, energy fiux 
above the 5 km s^^ wave amplitude decreases because a large fraction of the wave energy is 
dissipated in the chromosphere. When the wave amplitude becomes larger, the nonlinearity 
of the waves also becomes also larger. Then the waves experience frequent nonlinear mode 
coupling and dissipates rapidly. 

Figure [9] reveals the relationship between the transition region height and the amplitude 
of the waves. For the same reason described above, the spicule height is taller in the pink 
noise case than in the white noise case. 
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4.3. The waves that have the observed spectrum 



We have derived the photospheric horizontal velocity from the C-band movies by using 
LCT techniques. According to the figure [2|, the observed velocity spectrum generally has a 
double power law distribution. For the power spectrum density derived from the different 
regions listed in the Table [H we have performed a the double power law fitting. The fitting 
function is described as follows. 



The four parameters that define the spectral shape ( square root of the total power, 



Y< >, spectral break frequency, z/f,, and the two power index, aLiOtu) are described in 
the Table [H In order to distinguish the high and low frequency components, we implement 
the Fourier filter with a spectral break frequency. As a result, the low frequency components 
correspond to the granular motion while the high frequency components correspond to the 
inter granular motion. Since the power spectrum density in high frequency part is much 
smaller than that of low frequency part, waves driven by the granular motion is likely to be 
important for coronal heating. 

Then we have included the observed velocity spectrum listed in the first line of the 
Tabled] into our numerical simulations. In this case, the total energy fiux becomes 4.2 x 10^ 
erg cm~^ s~^ and the mean spicule height becomes 4.7 Mm. The results for the other case 
in the Table [T] are plotted in the figure [H] and IH] as white circles. Figure dO] is the Fourier 
spectra of the toroidal velocity with height. The color represents the power normalized by 
the total power at each height. There are at least three resonant peak at 1,3, and 5 mHz. 



Pu OC 




(18) 
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4.4. Synthesized spicules 

Although our numerical model is just 1 dimensional model along the flux tube, this 
flux tube has a three dimensional structure in real 3 dimensional space. Therefore we can 
synthesize the spicule motion in 3 dimensional space if we trace the fleld lines which have 
cool (less than 1 x 10^ K) plasma. In order to include projection effects, synthesized spicules 
are plotted in three different orientations with respect to an observer (see the right panel 
of flgure [TT]) . Although three bands of the spicules in flgure [TT] are synthesized from the 
same simulation, these spicules are viewed from different angles. The spicule with the white 
band is initially set on the plane that is made by the vertical direction and the line of sight 
direction, while the spicule with the thin shaded band is initially set on the sky plane. The 
spicule with the thick shaded band is initially set on the plane which is rotated 45 degree 
around vertical axis from the sky plane. 

We have compared the observed Ca II spicules with the synthesized spicules in the left 
panel of flgure [TTJ We superimposed the snapshot of the synthesized spicules on the Ca II 
image. Figure [11] is also available as an mpeg animation in the electronic edition of the 
Astrophysical Journal. From the movie, we can easily compare the height and the structure 
of spicules. The height of spicules from the photosphere deduced from the Ca II movie by 
eyes are comparable with that of synthesized ones. On the other hands, the movie clearly 
shows that the synthesized spicules are bent due to Alfven wave propagation, although the 
observed spicules usually have very straight structures. In the movie, white, thin shaded, 
and thick shaded spicules are changed into blue, black and red bands respectively. Although 
our model does not essentially depend on the foot point radius of the flux tube, appearance 
of the spicules can change. So we include the three different cases synthesized from three 
different foot point radius, rg = 15, 75, 150 km in the online movie. 
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5. Discussion 
5.1. Alfven wave resonant cavity 

Alfven waves are driven by the magneto-convection at the photosphere. The 



previous studies (1 


<!udoh & Shibata 


1999; 


Morivasu et al. 


2004; 


Morivasu & Shibata 


2004 


Suzuki & Inutsuka 


2005, 


2006 


) have assumed artificial wave drivers such as white noise 



spectrum because there were few observations of the photospheric temporal spectrum. In 
the present study, we include the observed disturbance into our numerical simulations. 

To investigate the dependence of coronal heating on the wave period, we have 
implemented a parameter survey on the wave period. As a result, we have found that the 
waves of period around 100-500 sec can transport the large amount of wave energy to the 
corona. These special period waves are also found to make the height of the transition 
region taller. Since the observed spectrum has large power around 300 sec, it is more 



efficie nt for coronal heating and spicule formation over white noise spectrum. 



Antolin et al. 



( 120081 ) also reported that the coronal temperature becomes highest when they generate the 



monochromatic Alfven waves of the period of 300 sec. 

The reason for the response of our model atmosphere against the special wave periods 
is considered to be Alfven wave resonance between t he ph otosphere and the transition 



region. Linear analysis by 



HollweeJ fll98lh and 



LeroyI (Il98ll ) suggested that the Alfven wave 



resonance can occur even when the fiux tube is not closed but opened. In the present 
study, our bottom boundary works as a refiective boundary since we have imposed the rigid 
boundary at the photosphere. The transition region also becomes a refiective boundary 
because the Alfven speed increases rapidly. Therefore standing waves appear in the region 
between the photosphere and the transition region and the standing waves will acquire 
larger amplitude by Alfven wave resonance. 
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Let us examine the linear properties of Alfven resonance in a flux tube in our model 
atmosphere by assuming a small wave amplitude using the same 1.5D MHD simulations. 
We assume small amplitude waves with white noise spectrum to check the Alfven wave 
resonance. Figure [12] is Fourier spectrum of the toroidal velocity and it clearly shows 
resonant behavior. The color represents the normalized power of the toroidal velocity at 
each height. The fundamental mode appeared around 1 mHz while the second harmonics 
appeared around 3 mHz. There are a lot of higher harmonics below 300 sec. Since the 
height of the transition region changes according to the wave property, resonant frequency 
will also change depending on the situation in the real Sun. When we include the observed 
velocity spectrum, the fundamental mode, ~ 1000 sec, has biggest power through the whole 
atmosphere. If the Alfven wave resonant cavity actually exists, the resonant periods are 
expected to be observed as spicule motion or coronal transverse velocity. 

In our simulation, the photospheric boundary is a fixed boundary by definition while 
the transition region works as a free boundary because of the high Alfven velocity at 
the corona. This means that the resonance occurs between the node (photosphere) and 
antinode (transition region). If, for simplicity, we assume the constant Alfven velocity (Va) 
in the resonance region , the resonance frequency Un can be described as follows. 



where L is the transition region height. When L equals to 2,000 km, and Va equals to 8 km 
s~^, resonance frequency becomes, 1,3,5 ... mHz. These frequencies agree quite well with 
our numerical experiments (Figure [T2|) . 

In the present study, we assume a reflective boundary at the photosphere. This 
assumption may not be valid if the flux tube extends straight down below the photosphere. 
However in the sub-photosphere, Alfven speed may be smaller than in the upper atmosphere 
and the flux tube can be easily bent by convection motion. In that case, the wave 




(19) 
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reflection can eventually occur at some point along the flux tube due to inhomogeneity 
of phase speed, resulting in a change in the resonant frequency. In order to confirm the 
boundary effect, we performed another simulation including the subphotospheric layer 
where the Alfven velocity decreases exponentially with depth. It was found that even if 
the additional subphotospheric layer exists, Alfven wave dynamics and energetics do not 
change qualitatively or quantitatively. 



5.2. Photospheric horizontal velocity spectrum 



We have investigated the granular motion by using LCT techniques for 14 different 
regions. LCT techniques require a stable data set with no effects from atmospheric seeing. 
The spatial resolution of Hinode/SOT is the highest of space optical observation of the 



Sun. For example. Solar Optic al Universa. 



telescope with 30 cm aperture (ITitle et al. 



Polarimeter (SOUP) aboard Spacelab 2 was a 



19861 ) while Hinode/SOT is a telescope with 50 



cm aperture. Although some ground-based observations have higher spatial resolution than 
that of Hinode/SOT, there are few data sets of this quality even if adaptive optics are used. 



of 



In general, t 



Tarbell et al. 



le re sults obtained in the present study are almost the same as that 



( 119901 ). We want to note that the derived LCT velocity is the apparent 
velocity, even though we have carefully removed the effect of p-mode. G-band intensity 
contrasts are considered to change due to the de struction of CH mo lecules when the plasma 



density decreases or the temperature increases (ISteiner et al. 



200l|). It is possible that 



the LCT velocity includes thes e intensity contrasts as well as the granular m otion. Since 



recent realistic 3D simulation (IStein fc Nordlund 



1998 



Carlsson et al 



2004 ) succeeded 



in reproducing granular motion, detailed comparison of LCT velocity with the simulated 
velocity should be investigated in the future. In order to avoid the ambiguity of L CT 



velocity, it is possible to calculate velocity by feature tracking (IBerger et al 



1998 
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MuUer et al.lll994[ ) but we do not consider this method. 



Another possibihties to generate the Alfven waves are rapid changes in the magnetic 
topology caused by reconnection or interac tion between flux tubes and acoustic waves com ing 



from ambient non-magnetic atmosphere (IHasan et aL 



2OO5I : 



Hasan fc van Ballegooijen 



20081 ). Recent observatio ns reveals that 



emerging in quiet region ( llshikawa et al. 



granu la r-sized trari s ient h orizontal flelds are 



20081). 



Isobe et al. 



fl2008[ ) found from 2D and 



3D MHD simulation that the high frequency Alfven waves can be created by magnetic 
reconnection between these transient flelds and ambient vertical flelds. 

Only the toroidal velocity is included in our simulations whereas the observed velocity 
is lateral. This may suggest that we should include the rotation of the observed velocity 
rather than lateral one. However, the basic equation of lateral motion is the same as that of 
toroidal motion except for the centrifugal force term. Since the centrifugal force turned out 
to be negligibly small in our simulations, it is valid to use troidal motion as a wave driver. 



5.3. 1.5 D approximation 



We have assumed a 1 dimensional model along the flux tube from the photosphere to the 
corona. Generally, there should be leakage of Alfven wave energy in the flux tube through 
the conversion to fast mode waves. Moreover, in the real Sun, there appears to be other 



mult i- dimensional eff ects such as Alfven fc fas t mode coupling by fleld bend ing (IWentze 



1974j ). phase mixing (IHeyvaerts fc Priestlll983l ). and Alfven wave resonance (jlonson 



19781). 



The synthesized spicules have bending structures due to Alfven wave whereas the observed 



spicules usually looks very straight. Indeed, wave lengths of kink Al: 



alon g spicules are e stimated to be more than 45,000 km ( iKim et al 



km ( iHe et al 



ven w aves propagating 



20081) or 2,000-4,000 



2OO9I ). Since the bending structures in our simulation are usually ~ 1,000 km. 
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there should be another dissipation mechanism of Alfven waves in multi-dimensional space. 



One possibility for multi-dimensional dissipation that can be suKgested by our 



1983h in the 



simulations is instabilities due to the standing waves (iHeyvaerts fc Priest 
Alfven wave resonant cavity. Shear flows at the antinodes of the standing waves can drive 
the Kelvin-Helmholtz instability, causing them to decay in a few periods. Moreover, current 
density at the nodes of the standing waves can cause the tearing mode instability. How 
much energy can be dissipated in such mechanisms will be interesting topics for future 
study. 

If there is greater Alfven waves dissipation in the chromosphere than the present 
model, the spicule heights must be reduced to some extent due to decreas e in nonlinearity . 



Although pure p-mode injection cannot explain the realistic spicule height (jSuematsu 



1990), 



p-mode may be c andidate for additional mo mentum in spicules. Slo w modes can 



3e created 



Takeuchi 



19991). 



by reconnectio n (ITakeuchi fc Shibata 



Nagata et al 



1998 



20011 ) or convective collapse (ISteiner et al. 
(l2008l ) observed the transient down flows with intensification 
of the field strength predicted by the model of convective collapse. These down flows are 
expect ed to be bounced back to the u pper atmosphere, resulting in a wave source for slow 



mode (jSterhng fc HoUweg 



1988 . 



19891). 



Recently, various simulations for nonlinear MHD wave propagation in mult i- dimensional 



space have been done. IShibatal (Il983l ) studied the nonlinear propagation of magnetoacoustic 
waves in the solar chromosphere in 1D/2D space and found that ID approximation is 
reasonable when the plasma beta is lower than unity. The propagatio n of nonlinear Alfv en 



waves were studied mechanism 



9851) . For coronal heating problem. 



or solar and astrop hvsic al iets (IShibata fc Uchida 



Hasan et al. 



(120051) and 



Hasan fc van Ballegooijen 



([20081) implemen ted the 2-D simulations of magnetoacoustic waves propagating in a flux 



tube. Moreover, IStein fc Nordlundl (119981 ) found the bending and horizontal displacement 
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of a flux sheet through the couphng with the convective motion self-consistently with 
radiative MHD code. These calculations only include the region below the chromosphere 
because the dynamical time scale becomes very short in the corona. Although recent 3D 
MHD simulations including the corona with fine spatial resolution can sufficiently resolve 
the wave dissipation and include the physics such as therm al conduction or radiative 



transfer (iGudiksen fc Nordlund 



2005 



Leenaarts et al. 



20091 ). it is worthwhile investigating 



ID simulation to understand nonlinear wave propagation. 



5.4. Energy equation 

When we want to precisely investigate the energy transfer between the photosphere and 
the corona, it is necessary to include radiative transfer and thermal conduction. The cooling 
time in the photo sphere and the chromosphere changes from 1 - 1,000 sec depending on 



height (IStix 



19701 ) while we assume a constant cooling time in the present study. Therefore 
the distribution of the temperature could be changed if we properly include radiative loss. 
Though the compressional mode waves are affected by radiative damping, it is inferred that 
non-compressional modes like Alfven waves are not significantly damped. Moreover, since 
the cooling time in the chromosphere (~ 200 sec) is much longer than the Alfven crossing 
time (~ 20 sec), the chromospheric dynamics of our simulation would not be greatly affected 
if radiative transfer was included. 

Thermal conduction is not so important as long as we are interested in the spicule 
dynamics. In our model, spicules are created by collision between transition region and 
slow mode shocks. The slow mode shocks are generated by nonlinear mode conversion from 
Alfven waves in the chromosphere. These formation process of spicules must not be affected 
by thermal conduction because the chromosphere is too cool. Thermal conduction become 
important when we want to study the energy balance in the corona or the transition region. 
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Suzuki fc Inutsukal (120051 . |2006| ) performed similar simulations to ours but including 
thermal conduction and radiative cooling. Their model succeeded in maintaining the 
corona and driving the fast solar wind by nonlinear low-frequency Alfven wave. Since the 
observed photospheri c velocity spectrurn is sirn ilar in shape and total power to the wave 



spectrum assumed in 



Suzuki fc Inutsukal (l2006l ). their model can maintain the corona when 



including the observed velocity spectrum. However the plasma beta in the photosphere 
and the chromosphere is very large in their simulation, therefore the dynamics, such 
as spicule motion, must be different. It is necessary to assume more realistic magnetic 
field distribution along the flux tube to explain spicule formation and coronal heating 
simultaneously. 



Carlsson &: SteinI (119971 ) derived the Ca II spectrum by using hydrodynamic simulations 
with non-LTE radiative transfer. In order to compare simulated results with observations, 
it is important to synthesize the spectrum during or after the simulations. Numerical 
simulations with radiative transfer are also needed when we solve the convection layer. 
Magnetic intensification through convective collapse and flux tube buffeting by convection 



motion can be included self-consistently by the radiative MHD simulations (iSteiner et al. 
19m . 



6. Conclusion 

We have investigated the nonlinear evolution of Alfven wave propagation along a flux 
tube in the solar atmosphere. In particular, we have concentrated on the dependence 
of the energy flux transmission and the spicule formation on the Alfven wave period. 
Monochromatic Alfven waves with given amplitude or Alfven waves that have continuum 
spectrum with various amplitudes are generated in our model atmosphere. Moreover we 
have analyzed the photospheric horizontal velocity spectrum from G-band movies and 
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included the spectrum in our numerical simulations. The results obtained in this paper are 
summarized as follows; 

1. The region between the photosphere and the transition region can be an Alfven wave 
resonant cavity. When including the observed velocity spectrum, we have obtained at 
least three resonant frequencies, 1,3, and 5 mHz. The Alfven waves that have resonant 
frequencies transfer more energy to the corona and increase the spicule height. 

2. Pink noise {Pi, oc u^^) has an advantage for coronal heating over white noise {Pi, oc u^). 
Alfven waves generated by pink noise also increase the spicule height. 

3. The photospheric horizontal velocity spectrum has larger power around 3 min than 
pink or white noise. The power around 3 min or more comes from granular motion 
while the high frequency components come from inter- granular motion. 

4. The observed velocity spectrum has an advantage for coronal heating over white 
noise. Obtained energy flux is 4.2 x 10^ erg cm~^ s~^, which means that Alfven waves 
transfer sufficient energy to heat the quiet corona. 

5. Structures of spicules synthesized from our numerical simulation are bending more 
than observed Ca II spicules. This implies the dissipation of Alfven waves by 
multi-dimensional effects that are not included in our simulations. 

Hinode is a Japanese mission developed and launched by ISAS/JAXA, with NAOJ as 
domestic partner and NASA and STFC (UK) as international partners. It is operated by 
these agencies in co-operation with ESA and NSC (Norway). Hinode is a Japanese mission 
developed and launched by ISAS/JAXA, collaborating with NAOJ as a domestic partner, 
NASA and STFC (UK) as international partners. Scientific operation of the Hinode mission 
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Fig. 1. — Example of G-band image observed by Hinode/SOT. The image in the upper 
right square (10"xlO") is zoom out image in the white square on the left. 
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Fig. 2. — Power spectrum density of the photospheric horizontal velocity. After the hor- 
izontal velocity is derived by the LCT technique, the spectrum is estimated by FFT. The 
solid line represents the mean power and dotted lines represent 1 sigma. 
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Fig. 3. — Structure of the flux tube that we assumed in our model. 
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Fig. 4.- 



Proflle of plasma beta with respect to height. 
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Fig. 5. — Dynamics of the Alfven waves in our typical simulation. Plots in each panel are 
stacked with time increasing upward in uniform increments of 8.0 s. Each panels represents 
the time profile of density (upper left), B<^/Bs (upper right), (lower left), and (lower 
right). 
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Fig. 6. — Total energy flux measured at the corona with respect to the wave period. The 
amplitude of each wave period is fixed to 1 km s~^. The white diamonds show each data 
point respectively while the vertical solid lines at each diamond represent the range of 1 
sigma. The horizontal dotted line indicates the sufficient energy to maintain the corona 
(3 X 10^ erg cm"^ s~-^). 
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Fig. 7. — Transition region height with respect to the wave period. The amphtude of each 
wave period is fixed to 1 km s^^. The white diamonds show each data point respectively 
while the vertical solid lines at each diamond represent the range of 1 sigma. 
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Fig. 8. — Total energy flux measured at the corona with respect to the wave amphtude. 
The sohd hne with white diamonds shows the white noise case while the dotted hne with 
asterisks shows the pink noise case. The white circle shows the observed spectrum case. The 
horizontal dashed hne indicates the sufficient energy to maintain the corona (3 x 10^ erg 
cm~^ s~^). 
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Fig. 9. — Transition region height with respect to the wave amphtude. The sohd hne with 
white diamonds shows the white noise case while the dotted hne with asterisks shows the 
pink noise case. The white circle shows the observed spectrum case. 
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Fig. 10. — Power spectrum density of toroidal velocity measured at each height. The power 
spectrum is normalized at each height so that the total power becomes one. The initial wave 
spectrum corresponds to the observed spectrum 7^1 hsted in Table [H 




Fig. 11. — Synthesized spicules compared with Ca II spicules observed by Hinode/SOT. 
Snapshot image from our simulation is shown in the upper right white square in the left 
panel. The distance between vertical white lines is 6 Mm. The right panel shows the initial 
positions of three synthesized spicules with respect to an observer. Three bands in the panel 
represent synthesized spicules that have different angle against the line of sight. 
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Fig. 12. — Power spectrum density of toroidal velocity measured at each height. The power 
spectrum is normalized at each height so that the total power becomes one. The Alfven 
wave is driven by white noise with 0.1 km s~^ amplitude. 
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Table 1: Horizontal velocity spectrum in quiet region. 



# of region 


Start Time 


Duration 


Mean Cadence 


\/< > 
y *r 




aL 


an 






(mill) 


(second) 


(km s-i) 


(niHz) 






1 


31-Oct-2006 11:10:31 


79 


30.0 


1.10 


4.7 


-0.6 


-2.4 


2 


02-NOV-2006 23:40:30 


99 


32.0 


1.07 


5.0 


-0.7 


-2.6 


3 


29-Dec-2006 16:50:30 


308 


30.0 


1.06 


4.3 


-0.6 


-2.1 


4 


30-Dec-2006 04:12:29 


272 


30.0 


1.07 


4.4 


-0.6 


-2.1 


5 


15-Feb-2007 08:14:50 


75 


34.2 


1.07 


4.8 


-0.6 


-2.4 


6 


17-Feb-2007 12:35:02 


174 


33.6 


1.11 


5.1 


-0.7 


-2.6 


7 


19-Feb-2007 18:19:02 


141 


33.6 


1.13 


5.0 


-0.7 


-2.5 


8 


27-Feb-2007 03:34:32 


119 


33.6 


1.07 


4.7 


-0.6 


-2.4 


9 


03-Mar-2007 07:35:31 


130 


33.6 


1.10 


4.8 


-0.6 


-2.5 


10 


09-Mar-2007 04:38:08 


81 


32.0 


1.10 


5.0 


-0.7 


-2.5 


11 


18-Mar-2007 07:56:04 


81 


32.0 


1.06 


5.0 


-0.7 


-2.6 


12 


06-Apr-2007 18:24:02 


345 


30.0 


1.08 


4.5 


-0.6 


-2.1 


13 


14-Apr-2007 18:13:31 


345 


30.0 


1.10 


4.5 


-0.6 


-2.1 


11 


29-Dee-2007 13:35:31 


278 


30.0 


1.13 


1.1 


-0.6 


-2.1 



